N93- 18838 


A COMPARISON OF VISCOELASTIC DAMPING MODELS 


Joseph C. Slater 
116 Jarvis Hall 
University at Buffalo 
Buffalo, New York 14260 

W. Keith Belvin 
Mail Stop 230 

NASA Langley Research Center 
Hampton, VA. 23665 

Daniel J. Inman 
1012 Furnas Hall 
University at Buffalo 
Buffalo, New York 14260 


Modem finite element methods (FEMs) enable the precise modeling of mass and 
stiffness properties in what were in the past overwhelmingly large and complex structures. 
These models allow the accurate determination of natural frequencies and mode shapes. 
However, adequate methods for modeling highly damped and highly frequency dependent 
structures did not exist until recently. The most commonly used method. Modal Strain 
Energy 1 * 2 , does not correctly predict complex mode shapes since it is based on the 
assumption that the mode shapes of a structure are real. Recently, many techniques have 
been developed which allow the modeling of frequency dependent damping properties of 
materials in a finite element compatible form. Two of these methods, the Gdla-Hughes- 
MeTavisht 3 " 4 method and the Lesieutre-Mingori 5 - 6 method, model the frequency dependent 
effects by adding coordinates to the existing system thus maintaining the linearity of the 
model. The third model, proposed by Bagley and Torvik 7 , is based on the Fractional 
Calculus method and requires fewer empirical parameters to model the frequency 
dependence at the expense of linearity of the governing equations. This work examines the 
Modal Strain Energy, Golla-Hughes-McTavish and Bagley and Torvik models and 
compares them to determine the plausibility of using them for modeling viscoelastic 
damping in large structures. 
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THE MODAL STRAIN ENERGY MODEL 


The most common method used for the modeling of viscoelastic damping in structures 
presently is the Modal Strain Energy method suggested by Ungar and Kerwin 2 . This 
method assumes that proportional damping {Rayleigh Damping) is an adequate model of 
the damping mechanisms of a structure. This implies that the modes of the damped 
structure are the same as that of the undamped structure. 

Modal Strain Energy begins with the complex stiffness representation of material 
damping properties. In this representation, the complex stiffness K* = fC + K"j where j 
represents the square root of - 1 , and fC and K" are the real and imaginary parts of the 
complex stiffness, respectively. The ratio K"/fC is the material loss factor. A more 
detailed description of the complex representation is given by Nashif, Jones and 
Henderson 8 . 

The loss factor of any mode i is given by the summation of the strain energy in each 
element, multiplied by its material loss factor, and divided by the total strain energy of the 
mode, i.e.. 



The variable if is the loss factor of the i th mode, rjj is the loss factor of the/* element 

at the i th natural frequency, V is the strain energy of the ft mode at a given amplitude, and 
Vj is the strain energy in the jth element when the structure is deformed in the ft mode 

shape at the same amplitude. The strain energy V in a structure or element with the 
stiffness matrix defined by K and the deformation defined by x is 


V = x r Kx 


( 2 ) 


Since the imaginary part of the global stiffness matrix is the assembly of the imaginary 
parts (AT') of the elemental stiffness matrices, equation (1) may be written 


rj = 


xf/Tx, 
x^fC x, 


( 3 ) 



where K is the real part of the global stiffness matrix and is denoted as K for the 
undamped and viscously damped systems. Note that this is precisely true in the case of the 
single degree of freedom system. This is a useful representation of the modal strain energy 
equation and will be referred to repeatedly. 

Although intuitively the concept of using energy ratios weighted by element loss 
factors is appealing, it has no theoretical basis. In the past there has not been an 
explanation of why modal strain energy is correct when the imaginary part of the stiffness 
matrix is proportional to the mass and stiffness matrices. It can be shown that modal strain 
energy is nothing more than the modal decoupling of a viscoelastic system where it is 
assumed that the imaginaxy part of the stiffness matrix must obey K'/w = C. 

The equations of motion for an unforced viscously damped multiple degree of freedom 
(MDOF) system may be written as 

My + Cy + Ky = 0 (4) 

Assuming a solution of the form 

y = u e imt (5) 

then substituting (5) into (4) gives 

-Mu?u + Ci(ou + Ku = 0 (6) 

The system equations written in terms of complex modulus corresponding to (4) are 

then 


My + (fC + K"i)y = 0 (7) 

Substituting (5) into (7) similarly gives 

-Mafu + (K' + K”i)u = 0 (8) 

Comparing (6) and (8), it is seen that for any given frequency, (Oj, 

C (Oj= K" 

which is the multiple degree of freedom representation of c=kr)/a>. 


( 9 ) 



likewise, substituting 


y = Pv (10) 

into (4), where P is the matrix of normalized eigenvectors of M ‘K and premultiplying by 
P r , equation (4) becomes 


PTMPv + P T CPv + P T KPv 

= diag ( m , )[ v + diag( 2 z { w, ) v + diag ( w f )v] - 0 (11) 

where are the modal masses. This is true if and only if C is proportional to M and K. 
Likewise, for the complex system, substituting (9) and (10) into (7), and premultiplying by 
/Ogives 

PTMP v + P T K' Pv + P T C (Oj P i v 


= diag(nti) [ v + ( diag(u), 2 ) + diag(2 & (oi) u)ji )v\- 0 (12) 

again if and only if C is proportional. This is identical to requiring that K' be 
proportional ( i.e., K"ko- aM+ pK). From ( 12) it can be seen that 

( P T FC P )-‘ ( P T C (OjP ) = diag(2 &) = diag(if) = H (13) 

when o)j = a However, using (9) gives 

( PTFC P )■> ( P T K" P ) = diag(2 &) = diag(rf) = H (14) 

Denoting the I th eigenvector as equation ( 14) may be written as 


which is identical to equation (3). Therefore, the same rules which apply to decoupling 
viscously damped systems apply to the proper use of modal strain energy. For non- 
proportionally damped systems, the matrix //defined by (14) will be non-diagonal and will 
give some indication of the non-proportionality of the system. 
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Thus, the modal strain energy technique is nothing more than the modal decoupling of 
a system with complex modulus damping. Also, the criteria used to define whether or not 
modal strain energy is a proper method for finding loss factors of a structure described by a 
complex modulus have been shown to be identical to those for decoupling a viscously 
damped system. 


THE GOLLA-HUGHES-MCTA VISH MODEL 


The Golla-Hughcs-McT avish (GHM) model is based upon the generalized standard 
linear model; however, it has been developed for direct incorporation into the finite element 
method. In the GHM model, the material complex modulus is written in the Laplace 
domain in the form 


E*(s) = E 0 (\ + h(s)) = 



(16) 


where the hatted terms arc free variables for curve fitting to complex modulus data and s is 
the Laplace domain operator. From ( 16) it can be seen that E*( to) - Eo for jco = 0 which 
means that no creep is allowed in this model. Also, the number of expansion terms, k, 
may be modified to represent the high or low frequency dependence of the complex 
modulus. In general between two and four terms are adequate. 


The finite element form of the GHM model for a single modulus and single expansion 
term is 


M 

0 

\(s) 

2 

0 

to 

z (s) 

s + 


0 

CO 



\\ + a)E 0 K -aE 0 R e A' 

’x(s) 


F(s) 

-aE^AX aE 0 A c 

z(s) 


0 


(17) 


where M is the original mass matrix, E n K is the original element stiffness matrix, Ac is 

the diagonal matrix of the non zero eigenvalues of K, and Rc is the matrix of the eigenvector 
associated with the eigenvalues of Ac- Even for the most complex linear elements, the 
GHM finite element remains linear and second order. Although the GHM system state 
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equations are much larger than the original undamped equations, the GHM state matrix is 
only slightly larger than the state matrix for a viscously damped system. Where the size of 
the state matrix for a viscously damped system is 2n x 2n, the state matrix of the GHM model 
is m x m where 


m= 



(18) 


Here n represents the number of viscoelastic elements, and pi and k\ represent the number 
of non-zero eigenvalues and number of expansion terms used in the /* element One 
drawback of this method which may be overcome simply is the addition of fictitious 
overdamped modes. These should be recognized as fictitious and discarded. 


FRACTIONAL CALCULUS - THE BAGLEY AND TORVIK MODEL 


The Bagley and Torvik fractional calculus viscoelastic model has been proposed based 
on the observations of Nutting^Gemant 10 ’ 1 1 , Caputo 1 -' 13 , Caputo and Minardi 14 , and 
Scott-Blair 15 that the mechanical properties of visooelastic materials seem to vary as a 
function of frequency raised to fractional powers. In the time domain, this represents 
fractional derivatives as defined by 


d° 

dt a 


x(t) = 


1 - f ^ dt 

T(1 - a) dtJ Q (t - x) a 


0<a< 1 


(19) 


where a represents the power of the derivative, F is the gamma function, and r is a 
dummy variable of integration. This in turn can be represented in the Laplace domain as 


£ 



= s a x(s ) 


( 20 ) 


where ^represents the Laplace transform operator. The general form of the fractional 
derivative model is then 


JL ,/A “ jPm 

a+ Jb n -?- r o=Ec+ 

dt Pn 


( 21 ) 
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The experimental results of Bagley and Torvik demonstrated that, for many materials, 
the stress-strain relation can be modeled well using oily the first expansion term in each 
series. In the Laplace domain, the Bagley and Torvik viscoelastic model is 

o(s) =- E °+ e(s) = fMs) e(s) (22) 

where fj(. s ) represents the complex modulus in the Laplace domain. In order to solve 
the final equations, a and /3 are restricted to fractional form. In the interest of brevity, no 
derivation is shown. One may be found in Bagley and Torvik 7 . 

The final form of the equations of motion are 

{B l s l/m +B 2 }y(s) = F(s) (23) 

where m is the smallest common denominator of a and fi, B\ and Ri are matrices of order 
nm(2 +b) x rvn ( 2 + b) , and y(s) and F(s) are the appropriate state vector and forcing 
function vector respectively. Equation (23) may be posed as an eigenvalue/eigenvector 
problem (setting F(s) =0) in order to solve for y (s) and s^ m . The system eigenvalues, s, 
and the system eigenvectors, y (s), may then be found using (23). However, notice that 
the order of the system is dramatically increased. For a second order viscously damped 
system, the size of the eigenvalue problem is 2n. The size of the Bagley and Torvik 
eigenvalue problem is [m(2+fl)] n. For the simplest possible Bagley and Torvik 
viscoelastic model, with a-b- 1/2, the order of the system is already Sn!. This will take 
6.25 times more memory and about four times as much time to calculate. Note that in the 
strictest sense this is not a finite element method, since no viscoelastic element has been 
developed which could be assembled into an existing finite element model in order to create 
a global FEM model. Another drawback of this method is the occurrence of unstable 
eigenvalues as described by Bagley and Torvik 7 . Although these may be disregarded when 
only interested in mode shapes and loss factors of modes, the forced response of this 
model would be unstable, which does not agree well with the real behavior of viscoelastic 
materials. 

It should be noted that much work has been done by Morgenthaler 16 of Martin 
Marietta using the Bagley and Torvik Model on the PACOSS program. The essence of his 
work is a numerical algorithm incorporating the accelerated subspace iteration technique to 
the complex modulus problem. Although the initial form of the stiffness matrix is assumed 
to be the fractional derivative, the algorithm’s first step is to evaluate the complex stiffness 



matrix at a frequency near the desired natural frequency. Then the desired natural 
frequency is found, the complex stiffness matrix is evaluated at the new frequency, and the 
desired natural frequency is found. The procedure iterates until the desired accuracy is 
reached, although it is mentioned that one iteration may be enough in many cases. The 
final state matrix is found by recoupling the eigenvalues and eigenvectors. 

This method discards the benefit of the Bagley and Torvik model by evaluating the 
frequency dependent stifiriess matrix instead of solving the complete set of frequency 
dependent equations derived by Bagley and Torvik 7 . This is not bad if all you are 
interested in are the correct mode shapes and eigenvalues. However, if that is the only 
goal, then there is no need to use the fractional derivative model. Any other Laplace 
domain representation which fits the modulus in the frequency domain would work equally 
as well. There is no reason not to simply use material data sheets directly to evaluate the 
complex stiffness matrix and avoid the curve fitting altogether. The end result of this 
method may model the system well, but it does not incorporate all of the frequency 
dependence of the materials in the final model. Any structural modification requires the 
complete recalculation of all of the desired eigenvalues and eigenvectors in order to find the 
new state space matrix, where the GHM method simply requires the assembly of a new 
element into the existing finite element model. 


AN EXAMPLE - THE EVOLUTIONARY VISCO-STRUT 


These three models have each been used to model a viscoelastic strut designed for use 
in the evolutionary model at NASA Langley. The Visco-Stmt is a load bearing member 
capable of supporting tensile and compressive forces in excess of 2000 lbs. In its present 
configuration, it has a static stiffness on the order of 30,000 lb/in. The viscoelastic material 
used is G.E SMRD, manufactured by the General Electric Astro Space Division. Both the 
GHM and the Bagley and Torvik models have proven to be capable of modeling the 
frequency dependent complex modulus of the Visco-Stmt well, while the Modal Strain 
Energy method simply uses the raw damping data for any element and therefore is not 
constrained by the need to curve fit The results of the curve fitting for the GHM and 
Bagley and Torvik models are shown below in figures 1 and 2. 
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For the GHM model, Kq = 3. 1 lx 10*. The remaining parameter are shown in Table 


1 . 



i=l 

i=2 

i=3 


2.29 

.520 

.319 

& 

9.83x10 14 

2.58x1023 

7.97x1020 

o>i 

1.91X10 18 

4.44x1025 

2.17X10 22 


Table 1. GHM Parameters for the Visco-Strut 



100 101 102 103 

Frequency (Radians/sec) 



Frequency (Radians/sec) 


Figure 1. Comparison of the GHM model (solid line) of the complex modulus and test 
data (dots). 



For the Bagley and Torvik model the parameters are A r (p2.86xl0 4 , ATj=3.6163xl0 3 , 
b= 1.9028x10-2 and a=p=l/2. 



10o 10i 102 103 


Frequency (Radians/sec) 

Figure 2. Comparison of the Bagley and Torvik model (9olid line) of the complex modulus 
and test data (dots). 

Some first attempts have been taken to model the effects of the Visco-Strut when 
placed in a small nine bay truss. To date, the only model which has been solvable is the 
MSE model. The Modal Strain Energy method is simple and quick because it uses the 
results from the dynamic model of the structure to find modal loss factors using equation 
(1) or (3). It has been shown by Johnson and Kienholz 1 to correctly find loss factors 
even when the damping is not proportional. Both the GHM and Bagley and Torv ik model 
solutions have encountered numerical difficulty. The GHM global FEM is ill-conditioned, 
while the size of the Bagley and Torvik model (1185 x 1 185) has caused significant 
numerical errors. Neither method has yielded useful results for this problem. 
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CONCLUSIONS 


The most robust method for determining modal loss factors is clearly the Modal Strain 
Energy method. If the system is proportionally damped (as determined by the matrix H in 
equation ( 14) being diagonal) then there is no need to use either the GHM or Bagley and 
Torvik method Even when the structure is non-proportional ly damped, there is little 
benefit to using either of the higher powered models unless there is a need to predictively 
model the mode shapes of the damped structure, or model the response of the structure to 
different excitations. However, if precise modeling of the response is necessary, one must 
decide whether the GHM modal will become too ill-conditioned for solution, or whether 
the large increase in the model size using Bagley and Torvik is acceptable. 
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